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TECHNICAL MEMORANDUM X- 


ROUGH-TO-SMOOTH TRANSIT lON OF AN 
EQUILIBRIUM NEUTRAL CONSTANT STRESS LAYER 


I. INTRODUCTION 

The purpose of the research reported herein is to develop a 
model of a low-level atmospheric flow over a terrain of changing rough- 
ness length. A flow of this kind is found at the windward end of landing 
strips as the wind blows onto the relatively smooth surface of the strip 
from the relatively rough surface of the surrounding terrain. The work 
involves use of the proposed model to derive functions which define the 
extent of the region and allow the calculation of wind and shear stress 
profiles at all points within the affected region. 

Elliott [l] was one of the first to study the change of terrain 
problem. He predicted that the thickness 6^ of the region affected by 

the change in roughness length increases with the 4/5-power of distance 
x from the discontinuity; i. e. , with fetch. Elliott refers to the region 
of the flow disturbed by the discontinuity as the internal boundary layer 
(IBL). Other investigators have pursued the problem using various 
models for the internal boundary layer; e.g., Panofsky and Townsend 
[2-5], Taylor [6], Peterson [7], and Rao, Wyngard, and Cote [8], 

Plate [9] has reviewed work done prior to 1971 on the internal boundary 
layer problem. The most recent studies utilize computer methods, 
together with a suitable turbulence closure model, to solve the govern- 
ing equations. 

A few investigators have studied the problem experimentally in 
wind tunnels. The earliest experimentation was done by Jacobs [ 1 0] , 
and the latest was accomplished by Antonia and Luxton [ 1 1 ] . The latest 
field data were reported by Bradley [ 1 2], Fair agreement has been 
obtained between experimental data and theoretical predictions. 

The technique of asymptotic matching in conjunction with assumed , 

scaling laws shall be used to extract as much information as possible 
from the equations of motions about the flow of concern. The motiva- j 

tion for this is to avoid lengthy and complicated numerical solution 
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methods which tend to mask the physics of the problem. The technique 
of asymptotic matching permits one to capture the physics of the prob- 
lem, extract basic needed information, and avoid the negative aspects 
of numerical methods. 


II. PROPOSED FLOW MODEL 

The present investigation deals with the transition region occur- 
ring in a neutral constant stress layer due to an abrupt change in surface 
roughness, the line of the discontinuity being normal to the direction of 
the wind. It is assumed that the surfaces are covered with uniformly 
distributed roughness elements (i. e., small obstacles or protrusions) 
the height of which may be characterized by the scale z^ (roughness 

length). The flow is assumed to be steady, two-dimensional, and 
incompressible, with neutral hydrostatic stability prevailing throughout 
the entire layer. The approach that is used in the analysis is similar 
to that of Clauser [13] in his study of the turbulent boundary layer and 
of Csanady [14] and Blackadar and Tennekes [ 1 5] in their investigations 
of turbulent Ekman layers. 

Blackadar and Tennekes [15] and many others have shown that 
the wind profile in the equilibri m surface layer upwind of the discon- 
tinuity is given by 


u 

u . 
*oi 



z 

z 

oi 


( 1 ) 


This equilibrium wind profile is represented schematically in Figure 1. 
Panofsky [16] states that the thickness of this surface, or logarithmic 
layer of the atmosphere, is of the order of 30 meters and that the 
Coriolis forces may be neglected in this layer. 

Previous investigators assumed that an internal boundary layer 

is formed above the smooth surface (roughness length of z ) and that it 

o 

grows with increasing distance from the discontinuity at x = 0. In 
Figure l two regions are defined by a curve emanating from the origin; 
i. e. , 


z = 6.(x) 


( 2 ) 
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Figure 1. Development of an internal boundary layer 
for transition of flow from rough to smooth terrain. 


which forms a surface above which the wind profile is described by 
equation (1) and below which the wind profile is undergoing transition to 
a new equilibrium profile given by 



(3) 


The fetch x corresponding to the achievement of this condition is repre- 
sented mathematically as x Equation (3) is of the same form as 

equation (1) and differs only in roughness length and velocity scale, 

z andu, , respectively, 
o vo 

A deficiency of the model of the flow as depicted in Figure 1 is 
that it does not provide for a gradual change of slope of the wind profile 
across the surface 6,(x). One can easily verify that the velocity gradient 

3u/0 z obtained from differentiating equation (1) is not equal to that 
obtained from differentiation of equation (3). Thus, although wind may 
be specified as continuous at the boundary of the two regions, the 
derivatives are not continuous across this boundary. 

To circumvent this difficulty it is proposed that the IBL be 
modeled as before but that an intermediate or matching region be 
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inserted between the two logarithmic regions as indicated in Figure 2. 


t 



Figure 2. Sublayers of the internal boundary layer 
of proposed model. 


The layer denoted as region I is the logarithmic layer in contact with 
the smooth surface. In region I the profile conforms to equation (3), 
with u a function of x. The thickness of this layer is denoted by 

with X < 6-« Region II is the new intermediate layer of thickness 6 - X 
i o 

in which a gradual transition from one logarithmic wind profile to the 

other is accomplished. A velocity defect la\% is hypothesized for 

region II. Along the surface defined by z = 5 , it is assumed that com- 

o 

plete matching to the external flow, governed by equation (1), has been 
accomplished, including wind, stress, and their derivatives. 

The thickness of the internal boundary layer as previously 

defined and as depicted in Figure 1 is intermediate between the outer 
layer defined by the surface z = 6 q (x) and the sublayer defined by the 

surface z = X(x). Along the surface z = 6.(x) the wind velocity u is 
given by 


4 


u. 

1 


( 4 ) 


u , k z 
*o 1 o 1 

This means that 6. defines the surface below which the flow has accele- 
i 

rated relative to the original equilibrium profile given by equation (1). 

Continuity requires that the region defined by 6. < z< 6 contain a flow 

1 o 

which has decelerated relative to the original equilibrium profile. 

Thus, the proposed model includes three length scales; viz., 

X, 6., and 6 , and three velocity scales; viz., u, , u,, and u, . 

i o " ! o: v v o 

These scales will be utilized in the next section to develop suitable 
similarity functions for wind and shear stress. 


III. PROPOSED SIMILARITY FUNCTIONS 

Townsend [3] found that similarity is possible in disturbed 
boundary layers if the disturbed layer is thin compared with the thick- 
ness of the boundary layer and if the roughness length is small compared 
with the thickness of the disturbed layer. Similarity may be assumed 
in the present problem if it is required that 6^ « h (where n is the 

thickness of the atmospheric boundary layer) and z q « 6^, which are 

clearly valid except verv close to and far from the origin. 

Two functions are defined: the velocity defect function F('n), 
which applies only to the intermediate layer (region II in Figure 2), and 
the shear stress function G(tj), which applies to the lower logarithmic 
layer (region I) as well as the intermediate layer. The velocity defect 
function is defined as 


u - u. 
F(q) = — 1 


(5) 


with the dimensionless height above the surface defined as 




z 


( 6 ) 
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The velocity u.is the velocity on the surface defined by z = j.(x) and is 
calculated from equation (4). Friction velocity u, in equation (5) is not 

the value at the surface but is a variable function of height z and fetch x. 

It is dependent on surface stress u* , a function of fetch only, and the 

>f o 

shear stress function G(r\), which is defined as 


G(ti) = 





(7) 


and which clearly indicates the relationship between u. v and G(t|). 

The nature of the similarity functions defined in equations (5) 
and (7) can be determined from a consideration of the boundary and 
matching conditions to which they must conform. The most obvious 
of these conditions are written first; i.e. , those for which the function 
is zero or unity. 

The function F(t)) is zero when the velocity u equals the velocity 
u^ on the surface defined by z = 6 Ax). The height of this surface is 

expressed nondimensionally as 



Thus, the condition at this surface can be written as 
F(/3) = 0 . 


(9) 


The function G(r|) is zero at the outer portion of region II where 
the shear stress must approach the upstream value , which is the 

value of shear stress found at all heights in the equilibrium layer 

upstream of the discontinuity. This condition occurs at the height 

z = 6 and is written as 
o 


G( 1 ) = 0 


( 10 ) 


? 

i. 


i 
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In the region above region lithe shear stress is constant; therefore, 

I 1 = 0 for z 2 6 (11) 

dz o 

The resulting condition for the shear stress functio. is 

G'(l) = 0 . (12) 

Similarly, in the lower part of the internal boundary layer; i.e. , in 
region I, 

G'fy) - 0 for t| — b , (13) 

where 



Conditions at the extremities of region I may be written as 


G'(b) = 0 


(15) 


and 


G ' (0) --- 0 . 


(16) 


In the region I the stress function is unity since t = t throughout the 

o 

layer. Thus 


G(r)| - 1 for t) 5 b , (17) 

Thus, all the boundary and matching conditions on the stress function 
are included in equations (10) through (17). Additional conditions on the 
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velocity defect function F(t]) will be sought next. 


At the outer limit of region II u = u , and equation (7) becomes 

o 


F(l) = - £ ln/3 


(18) 


where B is obtained from equation (6), u. from equation (4), and u from 

1 o 


u 

o 


u , 
:: oi 



(19) 


In order that the flow be self-preserving; i. e., that the function F(q) 
not depend on fetch x, it is required that F(l) = constant. This require- 
ment imposed on equation (18) implies that f3 is at most a function of 

z /z , which must be determined through experiments. The parameter 
o oi 

j3 will be determined later for a particular example. A condition on the 
derivative F'(l) may be obtained by matching derivatives at the juncture 
of region II and the upper logarithmic region. Equating expressions for 
3u/0z obtained from differentiation of equations (1) and (5) results in 
the condition 


F'(D = i 


( 20 ) 


Next the interfacial conditions as z — \ (r; -*• b) are considered. 
Again, the derivatives are matched, with the following result: 


F ' (b) = — 
1 kb 


( 21 ) 


Similarity requires that b, like B, is at most a function of z / z , 

o oi 

which must be determined by experiment. Its value is to be determined 
later for a particular case. 
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At this point a tentative velocity defect function is constructed 
which satisfies all of the necessary conditions for matching at the upper 
and lower bounds of region II; these conditions were given previously 
as equations (9), (18), (20), and (21). A suitable function for the region 
of velocity defect is 


F(t|) - £ In 



( 22 ) 


anti this function can be utilized to determine velocity distribution in 
region II after the growth law governing 6. and the constant fi are 

determined. The function given by equation (22) can be established via 
asymptotic considerations by demanding that the vertical gradient of u 
given by equation (3) as z -* °o and the vertical gradient of u given by 
equation (5) as z -* 0 merge together in an overlap region in which both 
formulations given by equations (3) and (5) are valid. This procedure 
assumes that z3u_,/3z becomes small in comparison to t)3F/3t| as 
z -* 0 


Equation (22) may be utilized to derive a relationship between 
surface friction and internal boundary layer thickness. At the outer 
interface of region I 


F(b) = 


u. 


u. 

i 


u , 
-■ : o 


(23) 


Substitution for u. from equation (4) and for u^ from equation (3) with 
z - A leads to 



u , 


where 



(25) 
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The utility of equation (24) will be demonstrated later. 

The required stress function must satisfy the matching condi- 
tions given by equations (10), (12), (15), (16), and (17). Figure 3 
depicts the shape of the required function as determined from the condi- 
tions listed above. An error function satisfies all matching conditions 
and is used to construct the stress function as is described below. The 
general form of the function is taken to be 


G(ti) = A erf (§) + C 


(26) 


where the argument £ is a function of q. 



G(tj) 

Figure 3. Schematic of properties of the 
nondimensional stress function G(t)). 


Let t| = /3 be chosen to be the point of inflection; i. e. , the maxi- 
mum absolute value of slope G'(r|). This choice follows from a con- 
sideration of the terms of the equation of motion which governs flow in 
region II; i. e. , 


3u 3u dr 
u -r— + w — = -r— 
ox o z d z 


(27) 
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It can be shown from equation (7) that the stress gradient is proportional 
to G , (‘ r l)« Equation (27) shows a direct dependence of stress gradient 
upon convective acceleration. Since the surface z = fi.(x) separates 

regions of net positive and negative acceleration, it is plausible that 
fluid in the vicinity of that surface experiences the largest accelerations 
and, hence, the largest stress gradients. Since maximum stress 
gradient implies maximum absolute value of G'(t)), then r\ = 0 corresponds 
to the maximum value of the derivative of the error function; i. e. , 

£ = 0 at t| = 0 . 

To construct the argument of £, cognizance is taken of the re- 
quirements stated in equations (12) and (15). The derivative of the 
error function is approximately zero when ||| > 2, where the error 
function itself is approximately unity. Choosing |£| = 4 for r| = b and 
T] = 1 and substituting these values in equation (26) enables a determina- 
tion of the constants C and A. To satisfy the conditions | = 4 at “n = b 
and | = 0 at t] = 0 , the argument 


P _ 4 (P--n) 

£ " 0 -b 


(28) 


is selected. Substituting £ = -4 at r| = 1 into equation (28) gives the 
relation between 0 and b; viz. , 


0 =| (b+1) 


(29) 


Substituting the above results into equation (26) yields the following 
expression for the shear stress function: 


G(q) = - erf 


' 4(1 b - 2 t|) 
1 -b 



(30) 


The constant b remains undetermined, and experimental data 
must be utilized to determine its value. Only one point of data is 
required, as will be demonstrated in the next section. 
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IV. DISTRIBUTION OF IBL THICKNESS AND 
SURFACE SHEAR STRESS 


In this section a relationship between the thickness 6. of the IBL 

and the fetch x, i. e. , the IBL growth law, is derived. This is obtained 
by integrating an asymptotic form of equation (27), the equation of 
motion. The terms are simplified by letting r| = fi and x -* 0. This 
results ir an easily integrated differential equation since F(/3), F'(J3) 
and u are eas ly inferred from equations (22) and (24) for the limiting 

values considei ;d. The surface friction velocity approaches zero in 
the limit as x -+ 0, which is consistent with equation (24). In view of 
this, it is plausible that u , at t] = j3 will approach zero as the discon- 
tinuity is approached; however, equation (7) indicates that a minimum 
value of u , is approached but not a zero value, because equation (7) is 

not valid at the roughness discontinuity. Nevertheless, in reducing the 
equation of motion it is assumed that u_, -*■ 0 as x -*• 0. A physical justi- 
fication for this assumption is that the velocity gradient, 


d u 
dz 


KJ '(£) u, ; 



(31) 


is expected to be small initially at r) = /3; i.e. , in the vicinity of the 
juncture of the new and old velocity profiles. This is expected because 
of the acceleration of the fluid inside the IBL relative to that above the 
IBL. Equation (11) requires a small value of u.,, as 6. shrinks near the 
origin. 


The above assumpt ^us reduce equation (27) to the simplified 
differential equation 

/«i \ 

6/ ln ^~ — j - -k z /3G'(/3) . (32) 


Integrati n of equation (32) results in the growth law. 


x 2 



r [ In r - (1 + M)] + e 


(33) 


M _ -k 2 /3G'(j3)x 
z 

o 


where 



o 


The boundary condition that r = e at x = o was used to derive 
equation (33). 

Equation (33) can be used to predict the growth of the internal 
boundary layer if the constant /3G'(,i) is determined from experimental 
data. From Bradley's data [12], valid for M = 4,83, as plotted in 
Figure 1 2 of Rao, Wyngaard, and Cote [8], 

j3G'(0) = - 1. 918 . (35 


It should be kept in mind that this numerical result is only valid for a 
particular value of z /z or rather M = 4. 83. 

o Ol 

The resulting value from equation (36) may be used with equation 
(33) to predict internal boundary layer thickness at other fetches. A 
curve depicting the results of this calculation is shown in Figure 4 with 
Bradley's data for M = 4. 83. 



Figure 4. Growth of the internal boundary layer with fetch. 
The experimental data are those of Bradley [ 1 2], 
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Values of 6^ calculated from equation (33) can be utilized with 

equation (24) to determine surface stress as a function of fetch. The 
curve in Figure 5 shows the comparison of predicted surface stresses 
with experimental data from Bradley [12] and numerical predictions 
from Rao, Wyngaard, and Cote [8], 



x (m) 

Figure 5. Nondimensional surface shear stress as a function of 
fetch. The experimental data are those of Bradley [12]. 

The numerical prediction of Rao, Wyngaard, and 
Cote [8] is indicated with the dashed curve. 

The prediction of the theory herein is 
indicated with the solid curve. 


The results shown in Figures 4 and 5 indicate the degree of 
success achieved in the prediction of internal boundary layer growth 
and of surface shear stress variation with fetch. Agreement between 
theoretical prediction and experiment is achieved at all station* for 
which experimental data are available. 
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V. PREDICTION OF WIND AND SHEAR STRESS 
PROFILES IN IBL 


The prediction of wind and shear stress profiles in the internal 
boundary layer following a discontinuity in surface roughness can be 
achieved through the use of the equations developed in previous sections. 
To achieve this, use is made of equations (22) and (30) in the form 


F(p) = 2. 5 In 



and 


GW . i 



(36) 


(37) 


in which the values b = 0.25 and /3 = 0. 625 have been substituted. The 
approximate value of b is obtained by differentiation of equation (26) and 
substitution of equation (35). The parameter /3 is evaluated with 
equation (29). 

To calculate a wind profile in the internal boundary layer it is 
necessary to know the two roughness lengths and the surface stress 
upwind of the discontinuity. With this information 6. is calculated from 

equation (33), u, from equation (24), and u. from equation (4). With 

'■'O x 

these values calculated as indicated above, a set of quantities is calcu- 
lated at each height z; for region II this procedure is carried out as 
follows: i) ri is calculated from equation (6); ii) G(r|) is calculated from 

equation (37); iii) u 2 (u 2 - t) is calculated from equation (7); iv) FOi) 

is (. alculated from equation (36); and v) u is calculated from equation 

(5). This completes the method of calculating wind and shear stress 

profiles in region II; i. e. , u - z and r - z profiles, respectively. 

Outside region II the wind profile is determined from equation (1) for 

z > g and from equation (3) for z < Some results of this calculation 
o 

are shown graphically in Figure 6 and are compared with the corre- 
sponding experimental data of Bradley [12], 
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Figure 6. Nondimens ional wind profile u + at x = 6. 1 m 
for M = 4. 83. The quantity u+ = u(z)/u(z = 112. 5 cm). 
The experimental data are those of Bradley [12]. 


It is evident from Figure 6 that the calculation method based on 
the proposed model is capable of producing reasonable predictions of 
wind profiles in the internal boundary layer associated with modification 
of the surface layer. Since the calculation of a wind profile requires 
use of the corresponding shear stress profile, the model is assumed to 
be valid for shear stress prediction as well. 
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